{-# LANGUAGE  TypeOperators, FlexibleInstances #-}
{- |
Module      :  Annealing
Description :  <optional short text displayed on contents page>
Copyright   :  (c) Sergey Sichevskiy 2012
License     :  BSD3

Maintainer  :  s.sichevskij@gmal.com
Stability   :  experimental
Portability :  portable

Простая реализация симуляции отжига методом Коши (стохастического метода
поиска глобального минимума функции в заданной области).

Метод отжига служит для поиска глобального минимума некоторой функции
f(x), заданной для x из некоторого пространства S, дискретного или
непрерывного.  Элементы множества S представляют собой состояния
воображаемой физической системы (энергетические уровни), а значение
функции f в этих точках используется как энергия системы E = f(x).  В
каждый момент предполагается заданной температура системы T, как правило,
уменьшающаяся с течением времени.  После попадания в состояние x при
температуре T, следующее состояние системы выбирается в соответствии с
заданным порождающим семейством вероятностных распределений G(x,T),
которое при фиксированных x и T задает случайный элемент G(x,T) со
значениями в пространстве S.  После генерации нового состояния x' = G(x,T),
система с вероятностью h(dE, T) переходит к следующему шагу в
состояние x0, в противном случае процесс генерации x0 повторяется.  Здесь
dE обозначает приращение функции энергии f(x')-f(x).  Величина h(dE,T)
называется вероятностью принятия нового состояния.

-}

module Annealing
  ( annealing
  )
    where

--import Text.Printf         (printf)

import System.Random       ( randomIO )
import Control.Applicative ( (<$>) )
import Data.HList hiding   ( hLength )

{--
The following pseudocode presents the simulated annealing heuristic. It
starts from a state s0 and continues to either a maximum of kmax steps or
until a state with an energy of emax or less is found.  In the process, the
call neighbour(s) should generate a randomly chosen neighbour of a given
state s; the call random() should return a random value in the range .  The
annealing schedule is defined by the call temperature(r), which should yield
the temperature to use, given the fraction r of the time budget that has
been expended so far.

  s ← s0; e ← E(s)                                  // Initial state, energy.
  sbest ← s; ebest ← e                              // Initial "best" solution
  k ← 0                                             // Energy evaluation count.
  while k < kmax and e > emax                       // While time left & not good enough:
    T ← temperature(k,kmax)                         // Temperature calculation.
    snew ← neighbour(s)                             // Pick some neighbour.
    enew ← E(snew)                                  // Compute its energy.
    if P(e, enew, T) > random() then                // Should we move to it?
      s ← snew; e ← enew                            // Yes, change state.
    if e < ebest then                               // Is this a new best?
      sbest ← snew; ebest ← enew                    // Save 'new neighbour' to 'best found'.
    k ← k + 1                                       // One more evaluation done
  return sbest                                      // Return the best solution found.
--}


{--
-- Функция Розенброка для двух переменных.
-- Она имеет глобальный минимум в точке (x,y) = (1,1) где f = 0.
f (x:::y:::Null) = (1-x)**2 + 100*(y-x**2)**2

m (x:::y:::Null) = (x':::y':::Null)
  where
    x' = 100 * exp (-x*x)
    y' = 100 - 100 * exp (-y*y)

main = do

 printf "x = %f\n" $ sqrt $ log (100::Double)
 annealing (f.m) (0:::0:::Null) 0.01 1500 1e-6
--}

--           Симуляция отжига методом Коши (быстрый отжиг).
--
-- Основным недостатком Больцмановского отжига является очень медленное
-- убывание температуры.  Например, для того, чтобы понизить исходную
-- температуру в 40 раз, требуется e^40 ~ 2.35 * 10^17 итераций, что уже
-- вряд ли приемлемо при решении каких-либо задач.  Ввиду этого Цу и Хартли
-- предложили алгоритм, который позволяет использовать для изменения
-- температуры схему
--
-- T(k) = T0/k

-- без потери гарантии нахождения глобального минимума. Это достигается за
-- счет использования в качестве G распределений Коши с плотностью
--
-- g(x', x, T) = T/(|x'-x|^2 + T^2)^((D+1)/2)
--
-- соответствующим образом нормированных. Например, в случае D = 1 приходим
-- к плотности
--
-- g(x', x, T) = 1/pi * T/(|x'-x|^2 + T^2).
--
-- К сожалению, это распределение не очень удобно моделировать в пространстве
-- размерности больше 1.  Этого можно избежать, например, с помощью
-- перемножения D одномерных распределений Коши:
--
-- g(x', x, T) = 1/pi^D * П T/(|x_i'-x_i|^2 + T^2), i = 1..D
--
-- но в этом случае нахождение глобального минимума гарантируется только при
-- законе изменения температуры не быстрее чем:
--
-- T(k) = T0/k^(1/D)
--
-- что гораздо медленнее схемы:
--
-- T(K) = T0/k.

annealing :: ( Show xs, HasLength xs, HasNeighbour xs ) => (xs -> Double) -> xs -> Double -> Int -> Double -> IO (xs,Double)
annealing f x0 t0 kmax emin = do

   (x',e') <- iteration 1 t0 (x0, f x0) (x0, f x0)

   return (x',e')

 where
  -- Выполнение одной итерации метода отжига
  iteration k t (x,e) (xb,eb) = do

     let t' = temperature k t

     (x' , e') <- oneStep t' (x,e)

     (xb',eb') <- if e' < eb 
                    then do
--                           printf "xy : %s\n" $ show (x',e')
                           return (x',e')
                    else do
                           return (xb,eb)

--     printf "iter : %d xy : %s e : %f t : %f\n" k (show x) e t'
--     printf "xy : %s\n" (show (xb',eb'))

     if k < kmax && e' > emin
        then do
               iteration (k+1) t' (x',e') (xb',eb')
        else do
               return (xb',eb')

  -- Переход в новое стояние
  oneStep t (x,e) = do

     x' <- neighbour x t

     let e' = f x'

     a  <- randomIO

     if moveProbability (e,e') t > a
        then do
               return (x',e')
        else do
               oneStep t (x,e)

  -- Принятая схема изменения температуры
  temperature i _ = t0 / k**n
    where
       k = fromIntegral i
       n = 1 / hLength x0

  -- Принятая вероятность принятия нового состояния
  moveProbability (e,e') t = 
    if (e'-e) < 0
      then 1
      else exp(-(e'-e)/t)


class HList xs => HasLength xs where
  hLength :: Num n => xs -> n

instance HasLength Null where
  hLength _ = 0

instance (HasLength xs) => HasLength (x:::xs) where
  hLength (_:::xs) = 1 + hLength xs



class HasNeighbour x where
  neighbour :: x -> Double -> IO x

instance HasNeighbour Null where
  neighbour _ _ = return Null

{--
instance (Random x, Floating x, Random y, Floating y) => C (x:::y:::Null) where
  neighbour (x:::xs) t = do

    a <- randomIO
    b <- randomIO

    let x' = tan $ pi * (a - 0.5)
        y' = tan $ pi * (b - 0.5)

    return $ x' ::: y' ::: Null
--}
  
instance (HasNeighbour xs) => HasNeighbour (Double:::xs) where
  neighbour (x:::xs) t = do

    a <- randomIO

    let x' = x + t * (tan $ pi * (a - 0.5))

    (x' :::) <$> (neighbour xs t)